• Go to github.com/jensroes/bristol-ws-2025
• Click on: Code > Download ZIP > unzip directory on your machine.
• Open project by double clicking on bristol-ws-2025.Rproj

• Provide a unified set of advanced statistical techniques that covers a wider range of data-analysis problems: multilevel generalised linear models
• General framework encapsulating a range of widely used techniques including linear regression, ANOVA, ANCOVA, MANCOVA, logistic regression, Poisson regression, log-linear, random effects / mixed-effects models etc.
• These models arise as special cases of a single underlying general principle.
• Apply this general principle to a wider range of data analysis problems.

• Linear regression: review of multiple linear regression and factorial ANOVA.
• General linear models: certain problems, which are poorly handled by ANOVA, can be easily accommodated (e.g. varying slope models).
• Generalised linear models: transformations applied to linear models allows modelling of data-types including categorical, ordinal, discrete frequency data.
• Multi-level generalised linear models: data with inherently hierarchical structure, which have often been inadequately dealt with.

• A model of how the probability distribution of one variable, known as the outcome variable, varies as a function of other variables, known as the explanatory or predictor variables.
• The most basic type of regression models is the normal linear model.
• In normal linear models, we assume that the outcome variable is normally distributed and that its mean varies linearly with changes in a set of predictor variables.
• By understanding the normal linear model thoroughly, we can see how it can be extended to deal with data and problems beyond those that it is designed for.

• We can use lm to easily implement such a model.

model <- lm(y ~ x, data = data)

• We can simulate data that exactly meet model assumptions (normal distribution, equal variance, independence etc.) with known model parameters.

# Set parameter values
n <- 1000
beta_0 <- 500 
beta_1 <- 50
sigma <- 100

# Random data for group 1
group_1 <- rnorm(n = n/2, mean = beta_0, sd = sigma)
# Random data for group 2
group_2 <- rnorm(n = n/2, mean = beta_0 + beta_1, sd = sigma)

# Generate data
sim_data <- tibble(group_1 = group_1,
                   group_2 = group_2) 

# Preview data
glimpse(sim_data)

Rows: 500
Columns: 2
$ group_1 <dbl> 558, 400, 470, 685, 448, 397, 456, 525, 530, 356, 440, 493, 49…
$ group_2 <dbl> 665, 562, 598, 459, 447, 606, 406, 564, 651, 393, 367, 639, 44…

# Change data format
sim_data <- pivot_longer(sim_data, cols = c(group_1, group_2), names_to = "x", values_to = "y")

# Preview data
glimpse(sim_data)

Rows: 1,000
Columns: 2
$ x <chr> "group_1", "group_2", "group_1", "group_2", "group_1", "group_2", "g…
$ y <dbl> 558, 665, 400, 562, 470, 598, 685, 459, 448, 447, 397, 606, 456, 406…

# As a reminder, these are the parameter values
beta_0 <- 500 
beta_1 <- 50
sigma <- 100

lm uses maximum likelihood estimation to determine for which set of parameter values are the data most likely.

# Normal linear model of simulated data
model <- lm(y ~ x, data = sim_data)

coef(model) # Model coefficients

(Intercept)    xgroup_2 
      498.6        43.4 

sigma(model) # Standard deviation

[1] 103

Check scripts

• exercises/normal_model_simulation_1.R
• exercises/normal_model_simulation_2.R

Observe how changing the number of observations affects the model estimates.

• In real life we don’t know the true parameter values.
• We also don’t know the true process that generates the data.
• For the simulated data we know the process that generated the data.
  • normal distribution because rnorm samples normally distributed data
  • both groups have the same variance because sigma was the same.
• We can use linear models to estimate parameter values from the data by which we make assumptions about the process that generates the data.

data <- read_csv("../data/martin-etal-2010-exp3a.csv") %>% glimpse()

Rows: 1,036
Columns: 4
$ item   <dbl> 11, 11, 11, 11, 11, 11, 25, 25, 25, 25, 25, 37, 37, 37, 37, 5, …
$ ppt    <dbl> 10, 9, 11, 6, 7, 2, 10, 9, 11, 7, 2, 10, 9, 6, 2, 10, 9, 11, 6,…
$ rt     <dbl> 1055, 2010, 461, 977, 1152, 1079, 1194, 1267, 1160, 684, 1367, …
$ nptype <chr> "conjoined", "conjoined", "conjoined", "conjoined", "conjoined"…

summarise(data, across(c(ppt, item), list(n = ~length(unique(.))), .names = "{col}"))

# A tibble: 1 × 2
    ppt  item
  <int> <int>
1    12    48

In standard repeated measures ANOVAs only one source of variance can be modelled at a time.

For example one would aggregate across items and by-participant neglecting the by-items variance.

data_ppt <- data %>% summarise(rt = mean(rt), .by = c(nptype, ppt)) 
data_ppt

# A tibble: 24 × 3
  nptype      ppt    rt
  <chr>     <dbl> <dbl>
1 conjoined    10  944.
2 conjoined     9 1298.
3 conjoined    11 1410.
4 conjoined     6 1000.
5 conjoined     7  922.
# ℹ 19 more rows

Each rt observation \(i \in 1 \dots N\) can be said to come a normal distribution with a mean \(\mu\) and a standard deviation \(\sigma^2\)

\[ \text{rt}_i \sim \mathcal{N}(\mu_i, \sigma^2)\\ \]

Using the linear regression equation, the mean \(\mu\) can be decomposed into

$$ _i = _0 + _1 _i + _i\

$$

where \(\text{nptype}_i\) is dummy coded (by default as treatment contrast in alphabetical order)

\[ \text{nptype}_i = \left\{ \begin{array}{ll} 0, \text{ if nptype}_i = \texttt{conjoined}\\ 1, \text{ if nptype}_i = \texttt{simple}\\ \end{array} \right. \]

Because

\(\mu_i = \beta_0 + \beta_1 \cdot 0 = \beta_0\),

\(\beta_0\) is the average rt for conjoined NPs.

Because

\(\mu_i = \beta_0 + \beta_\text{nptype} \cdot 1 = \beta_0 + \beta_1\) is the average rt for simple NPs, and \(\beta_1\) is the difference in rts between conjoined and simple NPs.

and finally \(\text{ppt}_i \sim \mathcal{N}(0, \sigma^2_\text{ppt})\)

library(afex)
m <- aov_car(rt ~ Error(ppt/nptype), data = data_ppt)
summary(m)

Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

              Sum Sq num Df Error SS den Df F value  Pr(>F)    
(Intercept) 28783750      1   587304     11  539.11 1.1e-10 ***
nptype          6217      1    13043     11    5.24   0.043 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

library(lme4)
m <- lmer(rt ~ nptype + (1|ppt), data = data_ppt)
anova(m)

Type III Analysis of Variance Table with Satterthwaite's method
       Sum Sq Mean Sq NumDF DenDF F value Pr(>F)  
nptype   6217    6217     1    11    5.24  0.043 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Complete exercise script linear-models/exercises/repeated_measures_model.R

and after that linear-models/exercises/mixed_effects_model.R

• Models with a combination of fixed (intercepts and slopes) and random effects.
• Repeated measures ANOVA are a special type of mixed-effects models.
• Linear models is not as contrained as ANOVAs.
• Fixed effects: systematic / deterministic effects (e.g. age, trial id, grouping variables)
• Random effects: non-systematic effects (e.g. ppts, items, class, country)
• Random effects can be nested: people within country
• Sometimes the differences can be blurry: comparing differences between two test sites rather than treating test sites as repeated measures variable.

• Repeated-measures ANOVA can only handle either participant or item variance
• Linear mixed-effects models can capture both (Baayen, Davidson, and Bates 2008)

m <- lmer(rt ~ nptype + (1|ppt) + (1|item), data = data)

summary(m)$coef

             Estimate Std. Error    df t value Pr(>|t|)
(Intercept)    1111.8       48.0  11.9    23.2 3.11e-11
nptypesimple    -33.4       14.5 977.9    -2.3 2.16e-02

• Likelihood ratio test

m_null <- lmer(rt ~ 1 + ( 1 | ppt ) + ( 1 | item ), data = data)
anova(m_null, m)

Data: data
Models:
m_null: rt ~ 1 + (1 | ppt) + (1 | item)
m: rt ~ nptype + (1 | ppt) + (1 | item)
       npar   AIC   BIC logLik deviance Chisq Df Pr(>Chisq)  
m_null    4 14319 14339  -7155    14311                      
m         5 14316 14340  -7153    14306  5.28  1      0.022 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Satterthwaite approximation for degrees of freedom (lmerTest)

m <- lmerTest::lmer(rt ~ nptype + (1|ppt) + (1|item), data = data)
summary(m)$coef

             Estimate Std. Error    df t value Pr(>|t|)
(Intercept)    1111.8       48.0  11.9    23.2 3.11e-11
nptypesimple    -33.4       14.5 977.9    -2.3 2.16e-02

• t-values larger than |2| roughly correspond to a significance level of \(\alpha = .05\) (Baayen 2008)
• Confidence intervals: 95% CI shows the range of hypotheses that can be ruled out for \(\alpha = .05\)

confint(m, "nptypesimple")

             2.5 % 97.5 %
nptypesimple -61.8  -4.93

• But who needs p-values anyway :)

summary(m)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: rt ~ nptype + (1 | ppt) + (1 | item)
   Data: data

REML criterion at convergence: 14289

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.518 -0.596 -0.140  0.448  4.965 

Random effects:
 Groups   Name        Variance Std.Dev.
 item     (Intercept)  1649     40.6   
 ppt      (Intercept) 25996    161.2   
 Residual             54536    233.5   
Number of obs: 1036, groups:  item, 48; ppt, 12

Fixed effects:
             Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)    1111.8       48.0   11.9    23.1  3.1e-11 ***
nptypesimple    -33.4       14.5  977.9    -2.3    0.022 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
nptypesimpl -0.151

summary(m)$coef

             Estimate Std. Error    df t value Pr(>|t|)
(Intercept)    1111.8       48.0  11.9    23.2 3.11e-11
nptypesimple    -33.4       14.5 977.9    -2.3 2.16e-02

(coefs <- m@beta) # extract estimates

[1] 1111.8  -33.4

# Predicted mean for np conjoined
coefs[1] + coefs[2] * 0

[1] 1112

# Predicted mean for np simple
coefs[1] + coefs[2] * 1

[1] 1078

# Check predictor contrasts (default: treatment)
contrasts(factor(data$nptype))

          simple
conjoined      0
simple         1

data_pred <- expand.grid(item = 1:48, ppt = 1:12, nptype = unique(data$nptype)) %>% 
  mutate(pred = predict(m, newdata = .),
         m = "random intercepts")

• Random slopes: NP-type effect is differently strong for certain participants and items.
• Correlation between random intercepts and slopes: E.g. NP-type effect might be weaker for participants / items with longer rts.
• Maximal random effects structure (Barr et al. 2013): for repeated measures within participants design we need random intercepts for participants and items with by-participants and by-items random slopes for NP type.

select(data, item, nptype) %>% unique() %>% arrange(item, nptype)

# A tibble: 96 × 2
   item nptype   
  <dbl> <chr>    
1     1 conjoined
2     1 simple   
3     2 conjoined
4     2 simple   
5     3 conjoined
# ℹ 91 more rows

select(data, ppt, nptype) %>% unique() %>% arrange(ppt, nptype)

# A tibble: 24 × 2
    ppt nptype   
  <dbl> <chr>    
1     1 conjoined
2     1 simple   
3     2 conjoined
4     2 simple   
5     3 conjoined
# ℹ 19 more rows

# Maximal random effects structure (with correlation)
m_max <- lmer(rt ~ nptype + (nptype|ppt) + (nptype|item), data = data)

# Maximal random effects structure (without correlation)
m_max_nocor <- lmer(rt ~ nptype + (nptype||ppt) + (nptype||item), data = data)

bind_rows(data_pred, data_pred_cor, data_pred_nocor) %>% 
  filter(ppt %in% c(7, 11)) %>% 
  mutate(across(m, ~case_when(str_detect(., "intercepts") ~ "\n(1|ppt) + (1|item)",
                              str_detect(., "without") ~ "\n(nptype||ppt) + (nptype||ppt)",
                              str_detect(., "with ") ~ "\n(nptype|ppt) + (nptype|ppt)"))) %>% 
  rename(`random effects` = m) %>% 
  ggplot(aes(x = nptype, y = pred, group = item)) +
  geom_point(size = .1) +
  geom_line(linewidth = .1) +
  facet_grid(ppt ~ `random effects`, scales = "free", labeller = label_both) +
  labs(y = "predicted rt", caption = "Lines represent different items")

# Maximal random effects structure (with correlation)
m_max <- lmer(rt ~ nptype + (nptype|ppt) + (nptype|item), data = data)
summary(m_max)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: rt ~ nptype + (nptype | ppt) + (nptype | item)
   Data: data

REML criterion at convergence: 14288

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.505 -0.596 -0.147  0.449  5.002 

Random effects:
 Groups   Name         Variance Std.Dev. Corr 
 item     (Intercept)   2749.8   52.44        
          nptypesimple   522.1   22.85   -1.00
 ppt      (Intercept)  27159.7  164.80        
          nptypesimple    52.4    7.24   -1.00
 Residual              54359.2  233.15        
Number of obs: 1036, groups:  item, 48; ppt, 12

Fixed effects:
             Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)    1111.9       49.3   11.5   22.58    7e-11 ***
nptypesimple    -33.6       15.0  125.7   -2.24    0.027 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
nptypesimpl -0.310
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Task: pair designs with their corresponding maximal random-effects structure!

• Design 1 (between participants / within items): both NP type conditions were presented for each item set but half of the participants saw simple NPs and the other half saw conjoined NPs.
• Design 2 (within participants / within items): each participant saw both conditions and each item was presented in both conditions.
• Design 3 (within participants / between items): half of all items were conjoined and the other half were simple NP type; each participant saw both condition (simple NPs and conjoined NPs).
• Design 4 (between participants / between items): half of all items were conjoined and half were simple NP type; half of the participants saw simple NPs and the other half saw conjoined NPs.

Model_A <- lmer(rt ~ nptype + (1|ppt) + (nptype|item), data = data) 
Model_B <- lmer(rt ~ nptype + (1|ppt) + (1|item), data = data)
Model_C <- lmer(rt ~ nptype + (nptype|ppt) + (nptype|item), data = data) 
Model_D <- lmer(rt ~ nptype + (nptype|ppt) + (1|item), data = data) 

Complete exercise script linear-models/exercises/mixed_effects_model_rirs.R

Type of outcome variable

• Continuous: lme4::lmer(dv ~ iv + (random slopes | random intercepts), data = data)
• Binary: lme4::glmer(dv ~ iv + (random slopes | random intercepts), data = data, family = binomial())
• Count: lme4::glmer(dv ~ iv + (random slopes | random intercepts), data = data, family = poisson())
  • Also negative-binomial lme4::glmer.nb and zero-inflated count models NBZIMM::glmm.zinb
• Ordinal: ordinal::clmm(dv ~ iv + (random slopes | random intercepts), data = data)

Frequentist mixed-effects models require familiarity with different packages.

Complete script linear-models/exercises/mixed_effects_model_brms.R

Frequentist mixed-effects models

library(lme4)
fit_lmer <- lmer(rt ~ nptype + (nptype|ppt) + (nptype|item), data = data)

Bayesian mixed-effects models

library(brms)
fit_brm <- brm(rt ~ nptype + (nptype|ppt) + (nptype|item), data = data)

• Two different ways of generalising from sample to the population:
• What do the data tell us about the population?
• Null-hypothesis significance testing
• Bayesian inference

• Classical / frequentist statistics; p-value based statistics
• p-value: How plausible are the data (or something more extreme) if we assume that there’s no effect?
• Evaluating the probability of data (e.g. t-value, F-statistic) assuming that the null hypothesis \(H_0\) is true.
• If the data are extreme enough, \(H_0\) is concluded to be implausible (rejected).
• If \(H_0\) is implausible, we assume \(H_1\) (alternative hypothesis).
• Statistically significant means data are unexpected under \(H_0\).

\[ Pr(\text{data} \mid H_0 = \text{TRUE}) \]

• Expressing uncertainty about a hypothesis expressed in probability distributions
• Updating one’s belief (prior, example later) about a hypothesis (e.g. difference between conditions, interactions etc) as the new information becomes available.

\[ Pr(H \mid \text{data}) \]

\[ \underbrace{Pr(\theta \mid \text{data})}_{\text{posterior}} \propto \underbrace{Pr(\theta)}_{\text{prior}} \cdot \underbrace{Pr(\text{data} \mid \theta)}_{\text{likelihood}} \] - Likelihood: probability of data given the model parameter value(s) - Prior: what did we already know about model parameter(s) - Posterior: our updated belief in the parameter value(s) after seeing the data

It tells us what we want to know!

• Summary stats of a Bayesian model are often very similar to what people think frequentist quantities (p-values, CIs) mean (Nicenboim and Vasishth 2016):
• What’s the relative evidence for one model (e.g., \(H_0\)) vs. another (e.g., \(H_1\))?
• What interval contains an unknown parameter value with .95 probability?

Instead of concluding

\[ p(\text{data} \mid H_0 = \text{TRUE}) \]

we can answer

\[ p(H \mid \text{data}) \]

• What’s the probability of an effect given the data (and what we already knew before)?
• Useful for small data sets for which prior information plays a huge role.
• Probabilistic sampling can handle convergence failures, missing data problems, complex models (many parameters, sources of random error), ceiling / floor effects, identifiability issues.
• Results are not dependent on sampling plan (Kruschke 2018).
• Flexible data modeling is now easily available (brms, rstanarm, rethinking).

• Estimate uncertainty over parameter values based on the observed data (Kruschke 2014).
• For this we need to start with

1. a probability model (e.g. a normal distribution and a set of predictors)
2. a prior (on e.g. means, difference, variability)

• brms uses probabilistic sampling to determine the posterior uncertainty about the parameter value.
• This involve calculating the weighted probability of the data in a potentially high dimensional parameter space (which can take time).

samples %>%
  filter(chain %in% 1, iteration %in% 1:100) %>%
  ggplot(aes(y=value, x = iteration, colour = chain)) +
  scale_color_viridis_d("Chains") +
  geom_path() +
  facet_wrap(~name, scales = "free") +
  labs(title = "Traceplot")

samples %>%
  filter(chain %in% 1:2, iteration %in% 1:100) %>%
  ggplot(aes(y=value, x = iteration, colour = chain)) +
  scale_color_viridis_d("Chains") +
  geom_path() +
  facet_wrap(~name, scales = "free") +
  labs(title = "Traceplot")

samples %>%
  filter(chain %in% 1:4, iteration %in% 1:100) %>%
  ggplot(aes(y=value, x = iteration, colour = chain)) +
  scale_color_viridis_d("Chains") +
  geom_path() +
  facet_wrap(~name, scales = "free") +
  labs(title = "Traceplot")

samples %>%
  filter(chain %in% 1:4, iteration %in% 1:250) %>%
  ggplot(aes(y=value, x = iteration, colour = chain)) +
  scale_color_viridis_d("Chains") +
  geom_path() +
  facet_wrap(~name, scales = "free") +
  labs(title = "Traceplot")

samples %>%
  filter(chain %in% 1:4, iteration %in% 1:500) %>%
  ggplot(aes(y=value, x = iteration, colour = chain)) +
  scale_color_viridis_d("Chains") +
  geom_path() +
  facet_wrap(~name, scales = "free") +
  labs(title = "Traceplot")

samples %>%
  filter(chain %in% 1:4, iteration %in% 1:2000) %>%
  ggplot(aes(y=value, x = iteration, colour = chain)) +
  scale_color_viridis_d("Chains") +
  geom_path() +
  facet_wrap(~name, scales = "free") +
  labs(title = "Traceplot")

samples %>%
  filter(chain %in% 1:4, iteration %in% 1000:2000) %>%
  ggplot(aes(y=value, x = iteration, colour = chain)) +
  scale_color_viridis_d("Chains") +
  geom_path() +
  facet_wrap(~name, scales = "free") +
  labs(title = "Traceplot")

• Traceplots: we want fat hairy caterpillars
• \(\hat{R}\) convergence statistic; should be \(<1.1\) (Gelman and Rubin 1992)
• Posterior predictive checks: compare observed data \(y\) and model predictions \(y_{rep}\)
• Complete script linear-models/exercises/model_diagnostics.R

• Prior knowledge about plausible parameter values.
• This knowledge is expressed as probability distributions (e.g. normal distributions).
• Help probabilistic sampler by limiting the parameter space.
• Small data samples are sensitive to prior information which makes intuitively sense.
• Otherwise data typically overcome the prior (automatic Ockham’s razor).
• Less common: test data against a (prior) effect suggested by the literature.

• brms uses default priors that can be changed as appropriate.
• Some defaults are good, (flat) priors are worth specifying.

Check defaults used earlier:

prior_summary(fit_brm)

p <- tibble(intercept = rstudent_t(10000, 3, 1039.5, 235),
       slope = runif(10000, -2000, 2000)) %>%
  ggplot(aes(x=intercept, y=slope)) +
    geom_point(size = .25, colour = "transparent") + 
  geom_density_2d_filled(alpha = 0.25, show.legend = F) +
  geom_density_2d(alpha = 0.5, show.legend = F, size = .15, color = "black") +
  scale_fill_viridis_d(direction = -1, begin = 0, end = .6) +
  scale_x_continuous(limits = c(0, 5000)) +
  scale_y_continuous(limits = c(-1000, 1000), breaks = seq(-1000, 1000, 250)) +
  labs(subtitle = "Parameter space")
ggExtra::ggMarginal(p, type="density", size=10, alpha = .25, margins = "both")

You will refit our model with more informative priors in the next script.

• Complete script linear-models/exercises/mixed_effects_model_brms_with_priors.R

• Summary of posterior parameter value estimates
• To test a hypothesis of parameter estimate (e.g. difference between groups, interactions, change in outcome):
  • ROPE
  • Savage-Dickey Bayes Factor: for comparison of nested models
• For comparison of models with different probability functions: Cross-validation (leave-one-out)

• Get posterior of slope \(\beta\) (difference between conditions)

beta <- as_draws_df(fit_brm) %>% pull(b_nptypesimple)

length(beta)

[1] 4000

beta[1:5]

[1] -34.0 -45.2 -33.0 -37.6 -23.7

• Determine ROPE to assess the uncertainty of effect size (Makowski, Ben-Shachar, and Lüdecke 2019)
• Calculate the standardised effect size \(\delta\) from posterior: \(\delta = \frac{\beta}{\sigma}\)
• ROPE of [-0.1, 0.1] is the range of negligible effect sizes (Cohn 1988; Kruschke 2018).

• Determine ROPE to assess the uncertainty of effect size (Makowski, Ben-Shachar, and Lüdecke 2019)
• Calculate the standardised effect size \(\delta\) from posterior: \(\delta = \frac{\beta}{\sigma}\)
• ROPE of [-0.1, 0.1] is the range of negligible effect sizes (Cohn 1988; Kruschke 2018).

sigma <- as_draws_df(fit_brm) %>% pull(sigma)
delta <- beta / sigma

abs(mean(delta)) # very small effect (see Cohen's d)

[1] 0.142

• Determine ROPE to assess the uncertainty of effect size (Makowski, Ben-Shachar, and Lüdecke 2019)
• Calculate the standardised effect size \(\delta\) from posterior: \(\delta = \frac{\beta}{\sigma}\)
• ROPE of [-0.1, 0.1] is the range of negligible effect sizes (Cohn 1988; Kruschke 2018).

rope(delta, range = c(-0.1, 0.1)) 

# Proportion of samples inside the ROPE [-0.10, 0.10]:

inside ROPE
-----------
26.34 %    

• The ROPE value indicates the extent to which the posterior cannot rule out a negligible effect.
• A meaningful effect size should have a small proportion of posterior samples within the ROPE.

Complete script linear-models/exercises/hypothesis_testing_rope.R

• Posterior: \(\mathcal{N}(-32.2, 17.4)\)
• Prior: \(\mathcal{N}(0, 100)\)

dnorm(0, 0, 100)

[1] 0.00399

dnorm(0, mean(beta), sd(beta))

[1] 0.00361

dnorm(0, 0, 100) / dnorm(0, mean(beta), sd(beta))

[1] 1.11

pd <- tibble(x = seq(-200, 200, by = 10),
             Prior = dnorm(x, 0, 100),
             Posterior = dnorm(x, mean(beta), sd(beta))) %>%
  pivot_longer(Prior:Posterior)

data_at_zero <- filter(pd, x == 0)
pr <- filter(data_at_zero, name == "Prior") %>% pull(value)
po <- filter(data_at_zero, name == "Posterior") %>% pull(value)

ggplot(pd, aes(x = x, y = value, colour = name)) +
  geom_line() +
#  geom_point(data = data_at_zero, aes(x = x, y = value), 
#             colour = "black", size = 2) +
#  geom_label(data = data_at_zero, aes(x = x, y = value, colour = name, label = round(value, 3)), 
#             colour = "black", size = 3, position = position_dodge(.0001)) +
  scale_color_colorblind("") +
  geom_vline(xintercept = 0, linetype = "dotted") +
  labs(y = "density", x = bquote(hat(beta))) #, subtitle = bquote("BF"==.(round(pr/po,2)))) +

• Posterior: \(\mathcal{N}(-32.2, 17.4)\)
• Prior: \(\mathcal{N}(0, 25)\)

dnorm(0, 0, 25)

[1] 0.016

dnorm(0, mean(beta), sd(beta))

[1] 0.00361

dnorm(0, 0, 25) / dnorm(0, mean(beta), sd(beta))

[1] 4.42

pd <- tibble(x = seq(-200, 200, by = 10),
             Prior = dnorm(x, 0, 25),
             Posterior = dnorm(x, mean(beta), sd(beta))) %>%
  pivot_longer(Prior:Posterior)

data_at_zero <- filter(pd, x == 0)
pr <- filter(data_at_zero, name == "Prior") %>% pull(value)
po <- filter(data_at_zero, name == "Posterior") %>% pull(value)

ggplot(pd, aes(x = x, y = value, colour = name)) +
  geom_line() +
#  geom_point(data = data_at_zero, aes(x = x, y = value), 
#             colour = "black", size = 2) +
#  geom_label(data = data_at_zero, aes(x = x, y = value, colour = name, label = round(value, 3)), 
#             colour = "black", size = 3, position = position_dodge(.0001)) +
  scale_color_colorblind("") +
  geom_vline(xintercept = 0, linetype = "dotted") +
  labs(y = "density", x = bquote(hat(beta))) #, subtitle = bquote("BF"==.(round(pr/po,2)))) +

• Posterior: \(\mathcal{N}(-32.2, 17.4)\)
• Prior: \(\mathcal{N}(0, 2000)\)

dnorm(0, 0, 2000)

[1] 0.000199

dnorm(0, mean(beta), sd(beta))

[1] 0.00361

dnorm(0, 0, 2000) / dnorm(0, mean(beta), sd(beta))

[1] 0.0553

pd <- tibble(x = seq(-200, 200, by = 10),
             Prior = dnorm(x, 0, 2000),
             Posterior = dnorm(x, mean(beta), sd(beta))) %>%
  pivot_longer(Prior:Posterior)

data_at_zero <- filter(pd, x == 0)
pr <- filter(data_at_zero, name == "Prior") %>% pull(value)
po <- filter(data_at_zero, name == "Posterior") %>% pull(value)

ggplot(pd, aes(x = x, y = value, colour = name)) +
  geom_line() +
#  geom_point(data = data_at_zero, aes(x = x, y = value), 
#             colour = "black", size = 2) +
#  geom_label(data = data_at_zero, aes(x = x, y = value, colour = name, label = round(value, 3)), 
#             colour = "black", size = 3, position = position_dodge(.0001)) +
  scale_color_colorblind("") +
  geom_vline(xintercept = 0, linetype = "dotted") +
  labs(y = "density", x = bquote(hat(beta))) #, subtitle = bquote("BF"==.(round(pr/po,2)))) +

• More informative priors give more weight to \(H_0\).
• Given such a prior, a posterior that favours \(H_1\) would be more convincing: BFs capture this.
• Regularising (skeptical) priors prevent the model to get overexcited by the sample (his words: McElreath 2020).
• If in doubt use weakly informative priors like student-t distributions that allow for larger tails, e.g. student_t(3, 0, 100)

• Explore hypothesis() in the exercise: linear-models/exercises/hypothesis_testing_bf_1.R
• Bonus: explore the Savage-Dickey method in linear-models/exercises/hypothesis_testing_bf_2.R

H: there is no difference between conjoined and simple NPs

hypothesis(fit_brm, "nptypesimple = 0")

Hypothesis Tests for class b:
          Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
1 (nptypesimple) = 0    -33.2      17.2    -66.3     0.75         NA        NA
  Star
1     
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

H: simple NPs are faster than conjoined NPs

hypothesis(fit_brm, "nptypesimple < 0")

Hypothesis Tests for class b:
          Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
1 (nptypesimple) < 0    -33.2      17.2    -61.4    -5.88         36      0.97
  Star
1    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

H: simple NPs are more than 10 msecs faster than conjoined NPs

hypothesis(fit_brm, "nptypesimple < -10")

Hypothesis Tests for class b:
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (nptypesimple)-(-10) < 0    -23.2      17.2    -51.4     4.12       11.1
  Post.Prob Star
1      0.92     
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

“Models are devices that connect theories to data. A model is an instantiation of a theory […]” (Rouder, Morey, and Wagenmakers 2016, 2)

• Our models describe how we understand reality.
• Model allows us to describe reality qualitatively (with their parameters) and quantitatively (parameter value estimates).
• Interpretation of parameters depends on data-modeling context.
• At minimum, distribution families (probability models) depend on data type.

• Gaussian: data come from normal distribution

fit <- brm(outcome ~ predictor + (1|participant), data = data, family = gaussian())

• Gaussian: data come from normal distribution
• Bernoulli: binomial outcomes (yes / no; correct / incorrect) (Winter and Bürkner 2021)

fit <- brm(outcome ~ predictor + (1|participant), data = data, family = bernoulli())

• Gaussian: data come from normal distribution
• Bernoulli: binomial outcomes (yes / no; correct / incorrect) (Winter and Bürkner 2021)
• Poisson: count data (number of words / mistakes)

fit <- brm(outcome ~ predictor + (1|participant), data = data, family = poisson())

• Gaussian: data come from normal distribution
• Bernoulli: binomial outcomes (yes / no; correct / incorrect)
• Poisson: count data (number of words / mistakes) (Winter and Bürkner 2021)
• Zero-inflated Poisson: count data with lots of zeros (number of typos per word)

fit <- brm(outcome ~ predictor + (1|participant), data = data, family = zero_inflated_poisson())

• Gaussian: data come from normal distribution
• Bernoulli: binomial outcomes (yes / no; correct / incorrect)
• Poisson: count data (number of words / mistakes) (Winter and Bürkner 2021)
• Zero-inflated Poisson: count data with lots of zeros (number of typos per word)
• Cumulative: ordinal (ordered categorical) data (Bürkner and Vuorre 2019) (Likert scales)

fit <- brm(outcome ~ predictor + (1|participant), data = data, family = cumulative()) # or acat, sratio

• Gaussian: data come from normal distribution
• Bernoulli: binomial outcomes (yes / no; correct / incorrect)
• Poisson: count data (number of words / mistakes) (Winter and Bürkner 2021)
• Zero-inflated Poisson: count data with lots of zeros (number of typos per word)
• Cumulative: ordinal (ordered categorical) data (Bürkner and Vuorre 2019) (Likert scales)
• Lognormal: zero bound data with positive skew; also skewed / shifted (log)-Normal, Wiener Diffusion models etc. [for RT data; Matzke and Wagenmakers (2009)]

fit <- brm(outcome ~ predictor + (1|participant), data = data, family = lognormal())

• Gaussian: data come from normal distribution
• Bernoulli: binomial outcomes (yes / no; correct / incorrect)
• Poisson: count data (number of words / mistakes) (Winter and Bürkner 2021)
• Zero-inflated Poisson: count data with lots of zeros (number of typos per word)
• Cumulative: ordinal (ordered categorical) data (Bürkner and Vuorre 2019) (Likert scales)
• Lognormal: zero bound continuous data with positive skew
• Ex-Gaussian: continuous data with positive skew

fit <- brm(outcome ~ predictor + (1|participant), data = data, family = exgaussian())

So far we used a Gaussian model which showed a poor fit to the data.

Alternative models: skewed normal, shifted (log)normal, Wiener Diffusion models, ex-Gaussian, mixture of lognormal etc. (see Matzke and Wagenmakers 2009; Roeser et al. 2024).

\[y \sim \mathcal{N}(\mu, \sigma^2)\]

• Data generating process can be described as a normal distribution with mean \(\mu\) and standard deviation \(\sigma\).

\[ \log(y) \sim \mathcal{N}(\mu, \sigma^2) \]

• Often used to address positive skew; e.g. response times (Baayen 2008)
• Log-values are only defined if \(y \in [0, \infty]\)
• De-emphasize large values over small values.
• Models the percentage change and not absolute differences:
  • Rather than saying “keystroke intervals are 40 msecs slower”, you say “keystrokes are 7% slower”.
  • Difference of 40 msecs can be huge for fast typing but negligible for pauses around 10 secs.

\[ \log(y - \exp(\text{ndt})) \sim \mathcal{N}(\mu, \sigma^2) \]

• Variation of the lognormal distribution: a zero shift would be equivalent to lognormal
• Includes a horizontal shift.
• Useful for modeling rt data because they are typically not just
  • right-skewed
  • non-negative
• But also
  • positively shifted: there is a minimum time required for any response

Accounts for a Minimum rt: The shift parameter represents the minimum possible rt, accounting for the fact that no response can be instantaneous.

Fit another probability model

• How accurately does the model predict new data.
• \(R^2\) et al.: more complex models are generally better ( \(R^2\) can’t go down).
• Overfitting: models with more parameters can make unreasonable predictions.
• Generalisations to data that were used to fit the model are over optimistic (Gelman, Hill, and Vehtari 2020).
• We will look at both fit and predictive performance of the model.

file_dest <- "../plots/fit2data.pdf"
include_graphics(file_dest)

• Cross validation: performance of a model on new data to remove the overfitting problem.
• Leave-One-Out (LOO) Information Criterion (LOO-IC):
  • Train model on \(N-1\) observations
  • Predict remaining data point from training model.
  • Repeat process \(N\) times to predict every observation from a model of the remaining data.
• Adding up prediction results gives an estimate of expected log-predictive density (\(elpd\)); i.e. approximation of results that would be expected for new data.
• loo() uses the probability calculations to approximate LOO-IC (i.e. Pareto smoothed importance sampling) (Vehtari, Gelman, and Gabry 2015, 2017).

loo(fit_brm)

Computed from 4000 by 1036 log-likelihood matrix.

         Estimate   SE
elpd_loo  -7141.4 40.8
p_loo        41.6  3.6
looic     14282.8 81.6
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.5]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

# product of n factors one for each data point: with theta being all parameters
# for LOO CV we perform exclude each data point one at a time which is equivalent to multiplying posterior by factor 1/p(y_i \mid \theta)
# LOO posterior excluding i is p(\theta \ mid y_{-i}) = p(\theta \mid y) / p(y_i \mid \theta) using the the 
# LOO distribution uses the posterior simulatios for \theta and giving each simulation a weight of 1 / p(y_i \mid \theta)
# Weighted simulations is used to approximate the predictice distribution of y_i (the held-out data point)
# $$p(\theta \mid y) \propto p(\theta) \prod^n_{i=1} p(y_i \mid \theta)$$

# pointwise density:
# Average likelihood of every observation in training sample
# this is done for each set of parameters sampled from the psoterior distribution
# Average liklihood for each obsersation
# Sum over all observations.
# resulting in the log-pointwise-predictive-density (lppd)

# effective number of parameters is the variance in log-likelihood for every observation: p_waic
# WAIC: -2*(lppd * p_waic)

log_lik(fit_brm) %>% 
  as_tibble() %>% mutate(sample = 1:n()) %>%
  pivot_longer(-sample, names_to = "obs", values_to = "loglik") %>%
  group_by(obs) %>%
  summarise(M = mean(loglik),
            var = var(loglik)) %>%
  mutate(obs = as.numeric(gsub(pattern = "V", replacement = "", obs))) %>%
  ggplot(aes(x = obs, y = M, ymin = M-var, ymax = M+var)) +
  geom_pointrange(fatten = .5) +
  labs(x = "Observations", y = "log-posterior predictive density (with variance)")

loo(fit_brm)

Computed from 4000 by 1036 log-likelihood matrix.

         Estimate   SE
elpd_loo  -7141.4 40.8
p_loo        41.6  3.6
looic     14282.8 81.6
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.5]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

#difference between two deviances has a chi-squared distribution; factor of 2 scales it that way; also called the Bayesian deviance


# N of parameters
# number of different conjectures for causes of explanations of the data
# How much iformation do we want the model to provide
# how flexible is themodel in fitting the training samples
# penalty term
# expected distance between in-sample and out-of-sample deviance

log_lik(fit_brm) %>% 
  as_tibble() %>% mutate(sample = 1:n()) %>%
  pivot_longer(-sample, names_to = "obs", values_to = "loglik") %>%
  group_by(obs) %>%
  summarise(M = mean(loglik),
            var = var(loglik)) %>%
  mutate(obs = as.numeric(gsub(pattern = "V", replacement = "", obs))) %>%
  ggplot(aes(x = obs, y = M, ymin = M-var, ymax = M+var)) +
  geom_pointrange(fatten = .5) +
  labs(x = "Observations", y = "log-posterior predictive density (with variance)")

#loo_R2(mixturemodel)
#bayes_R2(mixturemodel)

• Complete the script linear-models/exercises/model_comparison.R
• Which model has the lowest elpd score (i.e. highest predictive performance)?
• This script, again, assumes that you have stored the posterior of four models: Gaussian, lognormal, shifted lognormal, mixture of lognormals

mc <- readRDS("../models/model_comparison.rda")

mc$diff %>%
  as.data.frame() %>%
  rownames_to_column("model") %>% 
  select(model:se_elpd_loo) %>%
  mutate(across(where(is.numeric), round, 1)) %>%
  kable() %>%
  kable_styling("striped", full_width = F) %>%
  column_spec(1:5, width = "10em") %>%
  add_footnote("`elpd_loo` is reported as $\\\\widehat{elpd}$ and the difference as $\\\\Delta\\\\widehat{elpd}$.")

• brms isn’t more complicated than lme4 and comes with more flexibility.
• Priors require some thinking.
• Bayesian models allow straight forward interpretation of parameter values without an arbitrary threshold for significance.

• Bürkner tutorials: Bürkner (2018); Bürkner and Vuorre (2019); Bürkner (2019); Winter and Bürkner (2021)
• Vasishth tutorials: Sorensen, Hohenstein, and Vasishth (2016); Nicenboim and Vasishth (2016); Vasishth and Nicenboim (2016)
• Excellent book-length intros: Gelman, Hill, and Vehtari (2020), McElreath (2020), Kruschke (2014), Lambert (2018), Lee and Wagenmakers (2014)
• Convince yourself that Bayesian inference is the right answer: Clayton (2021)